A Continuous Parameter Homogeneous Semi-Markov Model for Stratigraphic Analyses from borehole Data

Markov Chain models are often used to represent subsurface stratigraphy and to simulate likely representations of the subsurface. The data used are cores or boreholes. Rock types are classified into a finite set of states, and a state transition probability matrix is estimated for some vertical sampling interval, ∆z. The estimated transition probability matrix is sensitive to the ∆z used. No objective methods for selecting a ∆z are available. The underlying deposition process is continuous rather than discrete. Hence, information is lost by discretizing the sampling domain, irrespective of the ∆z used. The thrust of this paper is to develop a Markov Chain model of subsurface stratigraphy that explicitly considers layer thickness statistics, recognizes the continuous nature of the underlying deposition process and does not suffer from the problems associated with sample discretization. A finite set of rock types is considered, and a state transition intensity matrix is estimated from interpreted borehole stratigraphic sequences. A simulation strategy is developed using the transition intensity matrix to determine the state by state transitions, and by bootstrapping (sampling with replacement from the empirical distribution function) a layer thickness corresponding to the new state from the borehole data. Example applications are provided to demonstrate the utility of the method. Continuous Parameter Homogeneous Semi-Markov Model by Ali and Lall 2 Introduction Stratigraphy in a sedimentary basin is usually characterized from well logs. Well log information is often classified into distinct rock or rock types. The vertical and horizontal stratigraphic changes are generally of subtle nature and can be probabilistically formulated and simulated in terms of such data. The goal of such a simulation is mainly the understanding and prediction of the stratigraphic behavior in a sedimentary basin (Harbaugh, and Bonham-Carter, 1970). The simulated stratigraphy helps locate stratigraphic traps of hydrocarbon, and search of oil, gas, or water in delineating reservoirs. Two types of probabilistic models for analyzing subsurface stratigraphy from such data are possible. First, one may consider an estimate of the unconditional probability of observing a particular rock type at a given location in the subsurface. Indicator Kriging (Journel, 1989) and Kernel intensity and regression estimators (Ali and Lall, 1993) have been used for this purpose. These methods work with a binary or categorical data set corresponding to rock types or hydraulic conductivity thresholds. Such estimates of the probability of occurrence of a rock type are useful for assessing the spatial variability of rock parameters over the site. However, they do not explicitly consider persistence in the occurrence of a particular type of rock, or the tendency of a certain rock type to follow another, e.g., coal seams commonly follow a seat-earth, and channel-conglomerates may be succeeded by point-bar sandstones (Miall, 1973). A Markov Chain (MC) approach is often used (Anderson and Goodman, 1957, Krumbein, 1968, Harbaugh, and Bonham-Carter, 1970, Bayer, 1985, and Sinvhal and Sinvhal, 1992) to assess the conditional probability of occurrence of different rock types through a one step state transition probability matrix. The utility of such an approach is to identify some genetic behaviors of the depositional process and generate synthetic stratigraphic profiles that possess such behaviors. Generating such profiles is useful for the understanding and prediction of the stratigraphic behavior in a sedimentary basin. Subsurface modelers may use such profiles for oil traps detection; and oil, gas, or water search in delineating reservoirs. Simulated stratigraphy also improves the understanding of flow behavior in stratified environments. Continuous Parameter Homogeneous Semi-Markov Model by Ali and Lall 3 Typically, MC models are applied for transitions along the depositional sequence, i.e, in the vertical. Transitions across states are usually considered across fixed stages or steps using a vertical sampling, ∆z, that is prescribed by the investigator. An appropriate ∆z depends on the scale of stratigraphic fluctuation, i.e., how rapidly rock types change in the vertical, and on the resolution of the borehole. A large ∆z relative to the scale of fluctuation leads to loss of information about the depositional sequence. On the other hand, if ∆z is taken to be too small, the number of transitions out of a given rock type or states may be very few, leading to a loss of the dependence structure and hence, to transition probability estimates that are not directly representative of the underlying process scale. This is a classical problem associated with discretization of an underlying continuous process. Despite much research (e.g., Sinvhal and Sinhval, 1992, Miall, 1973, Krumbein, 1967, Carr et al., 1966), a consensus for the optimal choice of ∆z in a given setting has not emerged. An alternative to such an approach is to use semi-Markov models (Krumbein and Dacey 1969, Dacey and Krumbein, 1970, Schwarzacher, 1972) where transitions between rock types occur at discrete time intervals, but the deposition rate of a particular rock type during such an interval is allowed to vary randomly. A continuous space model for lithology results. Past efforts at applying such models have suffered from a need to discretize the sample while estimating model parameters and/or assumptions of the parametric form of the probability distribution of rock thickness for each episode of deposition. In this paper we present a methodology for estimating parameters of a Homogeneous Markov Chain model in the vertical in a manner that obviates the need for a priori discretization of the domain. Transition intensities (relative frequency of state transitions per unit length) rather than transition probabilities are estimated from the data. A simulation strategy using these estimated transition intensities and a resampling of layer thicknesses is also developed. This approach is best classified as a semi-Markov Chain model. Applications to a data set from India and to one from Continuous Parameter Homogeneous Semi-Markov Model by Ali and Lall 4 Utah are presented to demonstrate the utility of the continuous parameter approach. Background A brief review of Discrete and Continuous Parameter Markov Chain models as used for stratigraphic modeling is offered in this section. Discrete Parameter Homogeneous Markov Chains (DHMC) Markov Chain (MC) models have been used for subsurface modelling since the 1950's. The occurrence of lithologies is viewed as a stochastic process. The lithology is modeled as a random state variable which takes discrete values X(z) (rock types: 1, 2, 3 ) as deposition, z, progresses. The values of the random variable X are called the states of the MC, while the deposition axis z is called the parameter of the MC. If the random variable X is considered to occur at discrete points (sampled at a resolution ∆z) along the z, axis, the Chain is called a discrete parameter MC. On the other hand, if the random variable X is allowed to occur at any point, the chain is called a continuous parameter MC. The Markovian property is stated as: P{X(zi) = k | X(zj) = l, X(z1) = m ...} = P{X(zi) = k | X(zj) = l} (z1< zj < zi) (1) where: zi, zj are the depths at points i, and j respectively. This property indicates that our knowledge of the state at point j is enough to infer the state at point i. First order dependence is expressed in equation (1). Models with higher order dependence can be considered. In reality, the deposition process varies continuously with time, and hence, with depth. The deposition type, rate, and its persistence vary over the site. However, a first order, homogeneous MC with a discrete parameter space (DHMC) has usually been used for the representation of such a process. The model is prescribed through a state transition probability matrix (TPM). The TPM is a two dimensional matrix which summarizes the relative transition frequencies Continuous Parameter Homogeneous Semi-Markov Model by Ali and Lall 5 from one state to another (Figure 1.). The TPM is computed from a tally matrix. Given n states, or types of rocks, an element fkm (∆z) in the tally matrix is the total number of transitions from state k to state m at a prescribed sampling interval ∆z. An element in the corresponding TPM, pkm(zi, zi+1), is the transition probability from state k at zi to state m at zi+1, and is given as: p ( z)=f z)/ f z) km km km m 1 n ∆ ∆ ∆ ( ( = ∑ (2) where ∆z is the sample discretization, (zi+1zi). A problem associated with the TPM calculation is to define an appropriate sampling discretization, ∆z. The sampling discretization, ∆z, is assumed to represent a lithologic unit within which the deposition is of the same type. There are two ways to discretize the sampling domain: 1) using break points for rock types, and 2) using a fixed ∆z. In the first approach, a sedimentological unit is presumed to have a uniform lithological composition, and transitions are considered across lithologic boundaries. Each unit, regardless of its thickness, represents a step for the MC. The lithologic thickness is not used. This approach may only be valid if the sedimentological units are of nearly the same thickness. The second approach provides an arbitrary increase in the number of observations. An equal interval sampling provides information about the lithologic thickness but it may result in information loss for lithologic transition. Bayer (1985, p91) showed that the tally matrix at a prescribed ∆z is equal to the sum of two matrices. The first matrix tallies the transitions to other states and has diagonal elements of value zero; while the second one tallies the lithologic frequency with off diagonal elements of value zero. The corresponding probability matrices represent transition (conditional), and total (unconditional) probability respectively. An illustration of TPM calculations for a synthetic situation with 3 types of rocks, with each layer of the same thickness (5 ft.), and a systematic alternation of layers is shown in Figure 1. For values of ∆z larger than or equal to the average bed thickness, (e.g., 5 ft. in Figure 1), there is a Continuous Parameter Homogeneous Semi-Markov Model by Ali and Lall 6 significant loss of information as many actual transitions are missed. For smaller values of ∆z, ( 2, or 1 ft. in Figure 1), the unconditional probability matrix becomes dominant in the TPM and the resulting transition probabilities become small. If the sampling interval is much smaller than the average bed thickness, the transition probability is extremely small (Davis, 1973, p. 285-286). If the beds have roughly the same thickness, a value of ∆z near the average thickness may provide a TPM that approximates the process to a certain extent. On the other hand, if the bed thickness is highly variable, a useful choice of ∆z may be difficult. This is usually the case in practice. This sampling problem was the focus of several studies. Krumbein (1967) suggested sampling intervals between 2 to 10 ft.. Sinvhal and Sinvhal (1979) suggested an interval of 2 m (6 ft.) or its multiple, and Miall (1973) found that a sampling interval slightly less than the average bed thickness gives satisfactory results. In general, no theoretical rules have emerged. Using the DHMC, simulation proceeds in a series of space steps ∆z based on the estimated TPM. Successive rock types are randomly sampled in steps of ∆z based on the TPM entries. There is usually no explicit consideration of layer thickness statistics. The thrust of this paper is to develop a Markov Chain model of subsurface stratigraphy that explicitly considers layer thickness statistics, recognizes the continuous nature of the underlying deposition process and does not suffer from the problems associated with sample discretization. Continuous Parameter Homogeneous Markov Chain and Semi-Markov Models Given the continuous nature of the deposition process, a continuous parameter homogeneous Markov Chain (CHMC) representation may be more appropriate than the DHMC. Keiding and Anderson (1989) define the transition intensity qkm as the ratio of the number of transitions from state k to state m to the total number of events of state k. Using this definition, the transition intensity qkm may be defined directly as the ratio of the total number of k-to-m transitions to the total thickness of state k observed along the entire profile of length L, and is given as: q n / L km km k * = (3) Continuous Parameter Homogeneous Semi-Markov Model by Ali and Lall 7 where nkm is the number of transitions from state k to state m, and Lk * is the total length of rock type k from which transitions to other rock types occur. Lk * does not include the last layer in each borehole from which no transition takes place. A Transition Intensity Matrix (TIM) comprised of elements qkm is used to describe a continuous parameter MC (Trivedi 1982, p. 362). A computation of the TIM using equation (3) for the synthetic data discussed earlier is shown in Figure 1(e). The quantity qkm is informative as it is a measure of both the lithologic transition and thickness. The direct definition of qkm in equation (3) is independent of ∆z. Krumbein (1968) used a continuous parameter Markov Chain to simulate stratigraphic sequences. The TPM is first estimated using a discretization ∆z. A TIM element qkm is then derived from the TPM as follows: